##### This script runs the linear mixed-effects models reported in the Table 2 and in the Table OA1.

### The required packages
library(lme4)
library(lmerTest)
library(performance)

##### The following set of models control for compositional effects (except for the null models - which do not include any explanatory variables)

### Model A: two levels with RE at the COUNTRY-ROUND level
# The null model
Model_A0 <-  lmer(red_sup ~ (1 | cntry_round), data = ess_4_ld_models, weights = pspwght)
summary(Model_A0)

# Model reported in table OA1
Model_A <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round), data = ess_4_ld_models, weights = pspwght)
summary(Model_A)
check_collinearity(Model_A)


### Model B: two levels with RE at the COUNTRY level
# The null model
ess_4_ld_models$cntry <- unclass(ess_4_ld_models$cntry)
Model_B0 <- lmer(red_sup ~ (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_B0)

# Model reported in table OA1
Model_B <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_B)
check_collinearity(Model_B)


### Model C: three levels with RE at ROUND level (3) and then COUNTRY-ROUND level (2)
# The null model
ess_4_ld_models$essround <- as.character(unclass(ess_4_ld_models$essround))
Model_C0 <- lmer(red_sup ~ (1 | cntry_round) + (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_C0)
icc(Model_C0, by_group = TRUE)

# Model reported in table OA1
Model_C <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_C)
check_collinearity(Model_C)


### Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
# The null model
Model_D0 <- lmer(red_sup ~ (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_D0)
icc(Model_D0, by_group = TRUE)

# Model reported in table OA1
Model_D <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_D)
check_collinearity(Model_D)


### Model E: cross-classified model
# The null model
Model_E0 <- lmer(red_sup ~ (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_E0)
icc(Model_E0, by_group = TRUE)

# Model reported in table OA1
Model_E <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_E)
check_collinearity(Model_E)


### Model F: full model
# The null model
Model_F0 <- lmer(red_sup ~ (1 | cntry_round) + (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_F0)
icc(Model_F0, by_group = TRUE)

# Model reported in table OA1
Model_F <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_F)
check_collinearity(Model_F)


### Model G: two levels with RE at the YEAR level
# The null model
Model_G0 <-  lmer(red_sup ~ (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_G0)
icc(Model_G0, by_group = TRUE)

# Model reported in table OA1
Model_G <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_G)
check_collinearity(Model_G)


# Plotting the model results
library(sjPlot)
tab_model(Model_A, Model_B, Model_C, Model_D, Model_E, Model_F, Model_G, digits = 3, digits.p = 4, digits.re = 4)
tab_model(Model_A0, Model_B0, Model_C0, Model_D0, Model_E0, Model_F0, Model_G0, digits = 3, digits.p = 4, digits.re = 4)



### Creating a table with ICCs for null models (i.e. models that only contain random intercepts and no explanatory variables)
# The following syntax creates Table 2 from the article
ICCs_null_models <- data.frame(c("Country", "ESS round", "Country-round"), matrix(data = rep(NA, 7*3), nrow = 3))
colnames(ICCs_null_models) <- c("Random effects", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# Model A: RE at cntry_round
ICCs_null_models[3, "Model A"] <- (as.data.frame(VarCorr(Model_A0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A0))[, "vcov"]))*100
icc(Model_A0, by_group = TRUE)

# Model B: RE at cntry
ICCs_null_models[1, "Model B"] <- (as.data.frame(VarCorr(Model_B0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B0))[, "vcov"]))*100
icc(Model_B0, by_group = TRUE)

# Model C: RE at essround and cntry_round
ICCs_null_models[2, "Model C"] <- (as.data.frame(VarCorr(Model_C0))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C0))[, "vcov"]))*100
ICCs_null_models[3, "Model C"] <- (as.data.frame(VarCorr(Model_C0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C0))[, "vcov"]))*100
icc(Model_C0, by_group = TRUE)

# Model D: RE at country and cntry_round
ICCs_null_models[1, "Model D"] <- (as.data.frame(VarCorr(Model_D0))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D0))[, "vcov"]))*100
ICCs_null_models[3, "Model D"] <- (as.data.frame(VarCorr(Model_D0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D0))[, "vcov"]))*100
icc(Model_D0, by_group = TRUE)

# Model E: RE at country and ess_round
ICCs_null_models[1, "Model E"] <- (as.data.frame(VarCorr(Model_E0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E0))[, "vcov"]))*100
ICCs_null_models[2, "Model E"] <- (as.data.frame(VarCorr(Model_E0))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E0))[, "vcov"]))*100
icc(Model_E0, by_group = TRUE)

# Model F: RE at all three effects
ICCs_null_models[1, "Model F"] <- (as.data.frame(VarCorr(Model_F0))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F0))[, "vcov"]))*100
ICCs_null_models[2, "Model F"] <- (as.data.frame(VarCorr(Model_F0))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F0))[, "vcov"]))*100
ICCs_null_models[3, "Model F"] <- (as.data.frame(VarCorr(Model_F0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F0))[, "vcov"]))*100
icc(Model_F0, by_group = TRUE)

# Model G: RE at essround
ICCs_null_models[2, "Model G"] <- (as.data.frame(VarCorr(Model_G0))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G0))[, "vcov"]))*100
icc(Model_G0, by_group = TRUE)

# Exporting Table 2 to a csv file
write.csv(x = ICCs_null_models, file = "table_2_ICCs.csv")

# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A0, Model_B0, Model_C0, Model_D0, Model_E0, Model_F0, Model_G0)



#################################################################################
# Storing and exporting results of the seven models into a Table OA1 of the paper #
#################################################################################

# Creating an empty data.frame into which the results will be stored
Models_initial <- data.frame(c("Intercept", "", "employed", "", "unemployed", "", 
                               "hinc_2_cop", "", "hinc_3_dif", "", "hinc_4_vdif", "", 
                               "female", "", "age", "", "isced2", "", "isced3", "", 
                               "isced4", "", "isced56", "", "stfeco", "", "lrscale", "", 
                               "LFPR_y0", "", "SOCEXP_y0", "", "GDP_y0", "", "gini_disp_y0", "", 
                               "liberal", "", "mediter", "", "Var_cntry", "Var_essround", 
                               "Var_cntry_round", "Var_individual", "ICC_cntry", "ICC_essround", "ICC_cntry_round"), matrix(data = rep(NA, 47*7), nrow = 47))
colnames(Models_initial) <- c("Parameter", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# loading the arm package (to be able to use the se.fixef function)
library(arm)
# Programming a function that will assign model coefficients (fixed effects) and respective standard errors to an article-ready table
assign_fe_results <- function(model_name, assignment_object, column_number){
  coefficients <- as.numeric(fixef(model_name))
  se_coefficients <- as.numeric(se.fixef(model_name))
  # Number of rows (for the variance components and other statistics) that will remain empty
  for(i in 1:((nrow(assignment_object)-7)/2)){
    assignment_object[(i*2)-1, column_number] <- coefficients[i]
    assignment_object[(i*2), column_number] <- se_coefficients[i]
  }
  return(assignment_object)
}

# Using the programmed function to create Table OA1 of the article
Models_initial <- assign_fe_results(model_name = Model_A, assignment_object = Models_initial, column_number = 2)
Models_initial <- assign_fe_results(model_name = Model_B, assignment_object = Models_initial, column_number = 3)
Models_initial <- assign_fe_results(model_name = Model_C, assignment_object = Models_initial, column_number = 4)
Models_initial <- assign_fe_results(model_name = Model_D, assignment_object = Models_initial, column_number = 5)
Models_initial <- assign_fe_results(model_name = Model_E, assignment_object = Models_initial, column_number = 6)
Models_initial <- assign_fe_results(model_name = Model_F, assignment_object = Models_initial, column_number = 7)
Models_initial <- assign_fe_results(model_name = Model_G, assignment_object = Models_initial, column_number = 8)
View(Models_initial)

### Manual assignment of variance components and ICCs
# Model A: RE at cntry_round
Models_initial[Models_initial$Parameter == "Var_cntry_round", "Model A"] <- as.data.frame(VarCorr(Model_A))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model A"] <- as.data.frame(VarCorr(Model_A))[, "vcov"][2]
icc(Model_A, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_cntry_round", "Model A"] <- (as.data.frame(VarCorr(Model_A))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A))[, "vcov"]))*100


# Model B: RE at cntry
Models_initial[Models_initial$Parameter == "Var_cntry", "Model B"] <- as.data.frame(VarCorr(Model_B))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model B"] <- as.data.frame(VarCorr(Model_B))[, "vcov"][2]
icc(Model_B, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_cntry", "Model B"] <- (as.data.frame(VarCorr(Model_B))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B))[, "vcov"]))*100

# Model C: RE at essround and cntry_round
Models_initial[Models_initial$Parameter == "Var_essround", "Model C"] <- as.data.frame(VarCorr(Model_C))[, "vcov"][2]
Models_initial[Models_initial$Parameter == "Var_cntry_round", "Model C"] <- as.data.frame(VarCorr(Model_C))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model C"] <- as.data.frame(VarCorr(Model_C))[, "vcov"][3]
icc(Model_C, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_essround", "Model C"] <- (as.data.frame(VarCorr(Model_C))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C))[, "vcov"]))*100
Models_initial[Models_initial$Parameter == "ICC_cntry_round", "Model C"] <- (as.data.frame(VarCorr(Model_C))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C))[, "vcov"]))*100

# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Models_initial[Models_initial$Parameter == "Var_cntry", "Model D"] <- as.data.frame(VarCorr(Model_D))[, "vcov"][2]
Models_initial[Models_initial$Parameter == "Var_cntry_round", "Model D"] <- as.data.frame(VarCorr(Model_D))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model D"] <- as.data.frame(VarCorr(Model_D))[, "vcov"][3]
icc(Model_D, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_cntry", "Model D"] <- (as.data.frame(VarCorr(Model_D))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D))[, "vcov"]))*100
Models_initial[Models_initial$Parameter == "ICC_cntry_round", "Model D"] <- (as.data.frame(VarCorr(Model_D))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D))[, "vcov"]))*100

# Model E: RE at country and ess_round
Models_initial[Models_initial$Parameter == "Var_cntry", "Model E"] <- as.data.frame(VarCorr(Model_E))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_essround", "Model E"] <- as.data.frame(VarCorr(Model_E))[, "vcov"][2]
Models_initial[Models_initial$Parameter == "Var_individual", "Model E"] <- as.data.frame(VarCorr(Model_E))[, "vcov"][3]
icc(Model_E, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_cntry", "Model E"] <- (as.data.frame(VarCorr(Model_E))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E))[, "vcov"]))*100
Models_initial[Models_initial$Parameter == "ICC_essround", "Model E"] <- (as.data.frame(VarCorr(Model_E))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E))[, "vcov"]))*100

# Model F: RE at all three effects
Models_initial[Models_initial$Parameter == "Var_cntry", "Model F"] <- as.data.frame(VarCorr(Model_F))[, "vcov"][2]
Models_initial[Models_initial$Parameter == "Var_essround", "Model F"] <- as.data.frame(VarCorr(Model_F))[, "vcov"][3]
Models_initial[Models_initial$Parameter == "Var_cntry_round", "Model F"] <- as.data.frame(VarCorr(Model_F))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model F"] <- as.data.frame(VarCorr(Model_F))[, "vcov"][4]
icc(Model_F, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_cntry", "Model F"] <- (as.data.frame(VarCorr(Model_F))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F))[, "vcov"]))*100
Models_initial[Models_initial$Parameter == "ICC_essround", "Model F"] <- (as.data.frame(VarCorr(Model_F))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F))[, "vcov"]))*100
Models_initial[Models_initial$Parameter == "ICC_cntry_round", "Model F"] <- (as.data.frame(VarCorr(Model_F))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F))[, "vcov"]))*100

# Model G: RE at essround
Models_initial[Models_initial$Parameter == "Var_essround", "Model G"] <- as.data.frame(VarCorr(Model_G))[, "vcov"][1]
Models_initial[Models_initial$Parameter == "Var_individual", "Model G"] <- as.data.frame(VarCorr(Model_G))[, "vcov"][2]
icc(Model_G, by_group = TRUE)
Models_initial[Models_initial$Parameter == "ICC_essround", "Model G"] <- (as.data.frame(VarCorr(Model_G))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G))[, "vcov"]))*100

# Exporting Table OA1 to a csv file
write.csv(x = Models_initial, file = "table_OA1_Models_composit_no_centering.csv")



### Collecting the coefficient estimates for the LFPR contextual variable and storing them to a data.frame
LFPR_initial <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                           matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_initial) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_initial[LFPR_initial$`Model type` == "Model A", 2:4] <- coef(summary(Model_A))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model B", 2:4] <- coef(summary(Model_B))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model C", 2:4] <- coef(summary(Model_C))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model D", 2:4] <- coef(summary(Model_D))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model E", 2:4] <- coef(summary(Model_E))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model F", 2:4] <- coef(summary(Model_F))[15, c("Estimate", "Std. Error", "df")]
LFPR_initial[LFPR_initial$`Model type` == "Model G", 2:4] <- coef(summary(Model_G))[15, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_initial$lower_bound <- LFPR_initial$Estimate - LFPR_initial$`Std. Error`*qt(p = 0.975, df = LFPR_initial$df)
LFPR_initial$upper_bound <- LFPR_initial$Estimate + LFPR_initial$`Std. Error`*qt(p = 0.975, df = LFPR_initial$df)



# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A, Model_B, Model_C, Model_D, Model_E, Model_F, Model_G)


